% Length Convergence Test

%% Calling Function to calculate V(Potential), R(Charge), Err(error of consecutive steps)
N = [501, 1001, 1501, 2001, 2501,4001];
L = 2000;
f = 1E-4;
% Physical Constant Setup (In SI Unit)
a0 = 0.5291772109217;  % Bohr Radius in Angstrom 
eps = 1;               % Relative permittivity
eps_0 = 8.8541878176E-12;   % Vacuum permittivity 
A = 2.62041966677E-20;      % Area per atom
vf = 1E6; 
e = 1.6E-19; % Electronic Charge
hbar = 1.05457172647E-34;
Vmax = sqrt(pi*f*(hbar^2)*(vf^2)/(A*e^2));
Rmax = f;
Vh = 27.2113850560;


%% Function Calling
for i = 1:length(N)
    [Va,Ra] = PNGraphene(f,N(i),L);
    
    % Data Conversion
    V = Va*Vh;
    R = Ra;
    
    % Data Arrangement
    N1 = length(V);
    x = linspace(-L*a0/2, L*a0/2,N1)'; 
    V = V';
    R = R';
    
    % Write to file
    filename = sprintf('pnGr(f=%1.1E,L=%i,N=%i)',f,L, N(i));
    dlmwrite([filename, 'V.txt'], [x,V]);
    dlmwrite([filename, 'R.txt'], [x,R]);
    
end

%% Data Plotting

N = [501, 1001, 1501, 2001, 2501, 4001];
L = 2000;
f = 1E-4;

figure; title('$\it{Basis Convergence Test}$');
subplot(2,1,1);
plot(-500:500,ones(1,1001)*Vmax, '--k'); hold on;
plot(-500:500,-ones(1,1001)*Vmax, '--k'); hold on;
xlabel('$\it{x(\AA)}$', 'interpreter', 'latex');
ylabel('$\it{V(x)(eV)}$', 'interpreter', 'latex');
subplot(2,1,2); hold on;
xlabel('$\it{x(\AA)}$', 'interpreter', 'latex');
ylabel('$\it{\sigma(x) (e/Atom)}$', 'interpreter', 'latex');

linespec={'-r', '-y', '-g', '-b', '-m','-k'};
for i = 1:length(N)
    filename = sprintf('pnGr(f=%1.1E,L=%i,N=%i)',f,L,N(i));
    display(filename);
    V = dlmread([filename, 'V.txt']);
    R = dlmread([filename, 'R.txt']);
    subplot(2,1,1); 
    plot(V(:,1),real(-V(:,2)),linespec{i});hold on;
    xlim([-10,10]);
    subplot(2,1,2); 
    plot(R(:,1),real(R(:,2)),linespec{i});hold on;
    xlim([-10,10]);
end
legend_str = strread(num2str(N,'N=%i\n'),'%s');
legend(legend_str);

    
    




